/*
  ivl_div.c  -  implementation of division functions

  This file is part of the IVL Library (LIBIVL).
  
  LIBIVL - library for interval computations 
  
  Copyright (c) 2007 Mateusz Paprocki

  The IVL Library is free software; you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation; either version 2.1 of the License, or
  (at your option) any later version.

  The IVL Library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with the IVL Library; see the file LICENCE.  If not, write to
  the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
  Boston, MA 02110-1301, USA.
*/

#include "interval.h"

int ivl_div(ivl_ptr ivl, ivl_srcptr x, ivl_srcptr y) {
	int inexact, inexact_left, inexact_right;
	
	inexact = inexact_left = inexact_right = 0;

	if (DOM_UNDEFINED(x) || DOM_UNDEFINED(y) || ivl_is_equal_u(y, 0))
		ivl->domain = IVL_UNDEFINED;
	else {
		switch (mpfr_sgn(IVL_LEFT(x))) {
			case SGN_NEGATIVE:
			
				switch (mpfr_sgn(IVL_RIGHT(x))) {
					case SGN_NEGATIVE:
					
						switch (mpfr_sgn(IVL_RIGHT(y))) {
							case SGN_POSITIVE:
							
								switch (mpfr_sgn(IVL_LEFT(y))) {
									case SGN_POSITIVE:
										inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), IVL_LEFT(y), GMP_RNDD);
										inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
										break;
									case SGN_ZERO:
										mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
										inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
										break;
									case SGN_NEGATIVE:							
										mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
										mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
										break;
								}
							
								break;
							case SGN_ZERO:
								inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
								mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								break;
							case SGN_NEGATIVE:
								inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
								inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), IVL_RIGHT(y), GMP_RNDU);
								break;
						}
					
						break;
					case SGN_ZERO:
						
						switch (mpfr_sgn(IVL_RIGHT(y))) {
							case SGN_POSITIVE:
							
								if (mpfr_sgn(IVL_LEFT(y)) > 0) {									
									inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), IVL_LEFT(y), GMP_RNDD);
									inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);
								} else {
									mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
									mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								}

								break;
							case SGN_ZERO:
								inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
								mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								break;
							case SGN_NEGATIVE:
								inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
								inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), IVL_RIGHT(y), GMP_RNDU);
								break;
						}
					
						break;
					case SGN_POSITIVE:
					
						switch (mpfr_sgn(IVL_RIGHT(y))) {
							case SGN_POSITIVE:
							
								if (mpfr_sgn(IVL_LEFT(y)) > 0) {									
									inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), IVL_LEFT(y), GMP_RNDD);
									inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDU);
								} else {
									mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
									mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								}

								break;
							case SGN_ZERO:
								mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
								mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								break;
							case SGN_NEGATIVE:
								inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDD);
								inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), IVL_RIGHT(y), GMP_RNDU);
								break;
						}
					
						break;
				}
			
				break;
			case SGN_ZERO:
			
				if (mpfr_sgn(IVL_RIGHT(x)) == 0) {									
					inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
					inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);
				} else {
					switch (mpfr_sgn(IVL_RIGHT(y))) {
						case SGN_POSITIVE:
						
							switch (mpfr_sgn(IVL_LEFT(y))) {
								case SGN_POSITIVE:
									inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
									inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDU);
									break;
								case SGN_ZERO:
									inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
									mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
									break;
								case SGN_NEGATIVE:
									mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
									mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
									break;
							}

							break;
						case SGN_ZERO:
							mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
							inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);
							break;
						case SGN_NEGATIVE:
							inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDD);
							inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);	
							break;
					}					
				}
			
				break;
			case SGN_POSITIVE:
			
				switch (mpfr_sgn(IVL_RIGHT(y))) {
					case SGN_POSITIVE:
					
						switch (mpfr_sgn(IVL_LEFT(y))) {
							case SGN_POSITIVE:
								inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), IVL_RIGHT(y), GMP_RNDD);
								inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDU);
								break;
							case SGN_ZERO:
								inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), IVL_RIGHT(y), GMP_RNDD);
								mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								break;
							case SGN_NEGATIVE:							
								mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
								mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
								break;
						}
					
						break;
					case SGN_ZERO:
						mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
						inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), IVL_LEFT(y), GMP_RNDU);
						break;
					case SGN_NEGATIVE:
						inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDD);
						inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), IVL_LEFT(y), GMP_RNDU);
						break;
				}

				break;
		}
				
		if (IVL_NAN_P(ivl))
			ivl->domain = IVL_UNDEFINED;
		else {
			if (ivl_has_zero(y))
				ivl->domain = IVL_PARTIAL;
			else
				ivl->domain = x->domain | y->domain;
			
			if (inexact_left)
				inexact += 1;
			if (inexact_right)
				inexact += 2;
		}
	}
	
	return inexact;
}

int ivl_div_u(ivl_ptr ivl, ivl_srcptr x, unsigned long y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || (y == 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		inexact_left = mpfr_div_ui(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
		inexact_right = mpfr_div_ui(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);		

		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_s(ivl_ptr ivl, ivl_srcptr x, long y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || (y == 0))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (y > 0) {
			inexact_left = mpfr_div_si(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_si(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
		} else {
			inexact_left = mpfr_div_si(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_si(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
		}

		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_d(ivl_ptr ivl, ivl_srcptr x, double y) {
	int inexact = 0, inexact_left, inexact_right;
	mpfr_t tmp;
	
	if (DOM_UNDEFINED(x) || (y == 0))
		ivl->domain = IVL_UNDEFINED;
	else {
		mpfr_init2(tmp, 53);
		
		if (y >= 0) {
			inexact_left = mpfr_set_d(tmp, y, GMP_RNDD);
			inexact_left |= mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), tmp, GMP_RNDD);
			
			inexact_right = mpfr_set_d(tmp, y, GMP_RNDU);
			inexact_right |= mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), tmp, GMP_RNDU);
		} else {
			inexact_left = mpfr_set_d(tmp, y, GMP_RNDD);
			inexact_left |= mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), tmp, GMP_RNDD);
			
			inexact_right = mpfr_set_d(tmp, y, GMP_RNDU);
			inexact_right |= mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), tmp, GMP_RNDU);
		}
		
		mpfr_clear(tmp);

		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_z(ivl_ptr ivl, ivl_srcptr x, mpz_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || (mpz_sgn(y) == 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		if (mpz_sgn(y) > 0) {
			inexact_left = mpfr_div_z(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_z(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
		} else {
			inexact_left = mpfr_div_z(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_z(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
		}
		
		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_q(ivl_ptr ivl, ivl_srcptr x, mpq_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || (mpq_sgn(y) == 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		if (mpq_sgn(y) > 0) {
			inexact_left = mpfr_div_q(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_q(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
		} else {
			inexact_left = mpfr_div_q(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
			inexact_right = mpfr_div_q(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
		}
		
		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_f(ivl_ptr ivl, ivl_srcptr x, mpfr_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || (mpfr_sgn(y) == 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		if (mpfr_sgn(y) > 0) {
			inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
			inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
		} else {
			inexact_left = mpfr_div(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
			inexact_right = mpfr_div(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
		}
		
		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_div_2exp(ivl_ptr ivl, ivl_srcptr x, unsigned long exp) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x))
		ivl->domain = IVL_UNDEFINED;
	else {		
		inexact_left = mpfr_div_2ui(IVL_LEFT(ivl), IVL_LEFT(x), exp, GMP_RNDD);
		inexact_right = mpfr_div_2ui(IVL_RIGHT(ivl), IVL_RIGHT(x), exp, GMP_RNDU);		

		IVL_CHECK(ivl, x->domain, inexact, inexact_left, inexact_right);				
	}
	
	return inexact;
}

int ivl_u_div(ivl_ptr ivl, unsigned long x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(y) || ivl_is_equal_u(y, 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		if (x > 0) {
			switch (mpfr_sgn(IVL_RIGHT(y))) {
				case SGN_POSITIVE:
				
					switch (mpfr_sgn(IVL_LEFT(y))) {
						case SGN_POSITIVE:
							inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
							inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
							break;
						case SGN_ZERO:
							inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
							mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
							break;
						case SGN_NEGATIVE:
							mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
							mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
							break;
					}
				
					break;
				case SGN_ZERO:
					mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
					inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
					break;
				case SGN_NEGATIVE:
					inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
					inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
					break;
					
			}
		} else {
			inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
			inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);
		}
		
		if (IVL_NAN_P(ivl))
			ivl->domain = IVL_UNDEFINED;
		else {
			if (ivl_has_zero(y))
				ivl->domain = IVL_PARTIAL;
			else
				ivl->domain = y->domain;
			
			if (inexact_left)
				inexact += 1;
			if (inexact_right)
				inexact += 2;
		}
	}
	
	return inexact;
}

int ivl_s_div(ivl_ptr ivl, long x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(y) || ivl_is_equal_u(y, 0))
		ivl->domain = IVL_UNDEFINED;
	else {		
		if (x > 0) {
			switch (mpfr_sgn(IVL_RIGHT(y))) {
				case SGN_POSITIVE:
				
					switch (mpfr_sgn(IVL_LEFT(y))) {
						case SGN_POSITIVE:
							inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
							inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
							break;
						case SGN_ZERO:
							inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
							mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
							break;
						case SGN_NEGATIVE:
							mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
							mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
							break;
					}
				
					break;
				case SGN_ZERO:
					mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
					inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
					break;
				case SGN_NEGATIVE:
					inexact_left = mpfr_ui_div(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
					inexact_right = mpfr_ui_div(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
					break;
					
			}
		} else if (x < 0) {
			switch (mpfr_sgn(IVL_RIGHT(y))) {
				case SGN_POSITIVE:
				
					switch (mpfr_sgn(IVL_LEFT(y))) {
						case SGN_POSITIVE:
							inexact_left = mpfr_si_div(IVL_LEFT(ivl), x, IVL_LEFT(y), GMP_RNDD);
							inexact_right = mpfr_si_div(IVL_RIGHT(ivl), x, IVL_RIGHT(y), GMP_RNDU);
							break;
						case SGN_ZERO:
							mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
							inexact_right = mpfr_si_div(IVL_RIGHT(ivl), x, IVL_RIGHT(y), GMP_RNDU);
							
							break;
						case SGN_NEGATIVE:
							mpfr_set_inf(IVL_LEFT(ivl), SGN_NEGATIVE);
							mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
							break;
					}
				
					break;
				case SGN_ZERO:
					inexact_left = mpfr_si_div(IVL_LEFT(ivl), x, IVL_LEFT(y), GMP_RNDD);
					mpfr_set_inf(IVL_RIGHT(ivl), SGN_POSITIVE);
					break;
				case SGN_NEGATIVE:
					inexact_left = mpfr_si_div(IVL_LEFT(ivl), x, IVL_LEFT(y), GMP_RNDD);
					inexact_right = mpfr_si_div(IVL_RIGHT(ivl), x, IVL_RIGHT(y), GMP_RNDU);
					break;
					
			}
		} else {
			inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
			inexact_right = mpfr_set_ui(IVL_RIGHT(ivl), 0, GMP_RNDU);
		}
		
		if (IVL_NAN_P(ivl))
			ivl->domain = IVL_UNDEFINED;
		else {
			if (ivl_has_zero(y))
				ivl->domain = IVL_PARTIAL;
			else
				ivl->domain = y->domain;
			
			if (inexact_left)
				inexact += 1;
			if (inexact_right)
				inexact += 2;
		}
	}
	
	return inexact;
}

